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We study the dynamics of perturbations in time delayed dynamical systems. Using a suitable 
space-time coordinate transformation, we find that the time evolution of the linearized perturbations 
(Lyapunov vector) can be described by the linear Zhang surface growth model [Y.-C. Zhang, J. 
Phys. France 51, 2129 (1990)], which is known to describe surface roughening driven by power-law 
distributed noise. As a consequence, Lyapunov vector dynamics is dominated by rare random events 
that lead to non-Gaussian fiuctuations and multiscaling properties. 
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Deterministic chaos is an ubiquitous feature of nonlin- 
ear systems and much is understood about its dynam- 
ical properties in low-dimensional systems. Lyapunov 
exponents, measuring the separation of initially nearby 
trajectories, have demonstrated to be helpful in charac- 
terizing chaotic motion in systems with a few degrees 
of freedom 0. In the case of spatially extended non- 
linear systems, a straightforward generalization leads to 
the concept of Lyapunov exponent density or spectrum, 
but the problem here is more complex and one must take 
into account propagation and diffusion phenomena that 
complicate much theproblem 0. 

Pikovsky et al. |j| have recently shown that the 
Lyapunov vector, which contains the spatial dependence 
of the diverging growth rate of the linearized perturba- 
tions, can be described as a nonequilibrium growing inter- 
face. Moreover, they have reported that in many cases, 
as for instance in coupled-map lattices, Lyapunov vector 
growth is consistent with that of the well-known Kardar- 
Parisi-Zhang (KPZ) equation of surface growth. In 
the case of systems with a conservation law, one has to 
be more cautious since KPZ scaling might be not satis- 
fied, e.g., for the Hamiltonian lattices studied in Ref. [^. 
An interesting question that naturally arises is whether 
the dynamics of perturbations in extended systems could 
be divided in a few universality classes of growth accord- 
ing to the existence of symmetries and/or conservation 
laws, as occurs in nonequilibrium surface growth 

In this Letter we study Delayed Dynamical Systems 
(DDSs), which are formally infinite dimensional dynami- 
cal systems and show many aspects of space-time chaos. 
DDSs appear in a number of biological and physical situ- 
ations, e.g., regulation of blood cells production lasers 
with delayed feedback 0,0, etc. In recent years, DDSs 
have attracted a renewed interest as they have been pro- 
posed as candidates for secure communication systems 
based upon semiconductor lasers with time-delayed opti- 
cal feedback [HEl- 

Our aim in this Letter is to show that the dynamics 
of the Lyapunov vector in DDSs can be generically de- 



scribed by a stochastic surface growth equation in which 
rare fluctuations play a most important role. We ar- 
gue that the growth model belon gs t o the universality 
class of the linear Zhang model [Jill [11 [11 13 . This 
is in contrast with other spatially extended dynamical 
systems that in general lead to KPZ-like growth of the 
amplitude perturbations [1 13| . Our results clearly show 
that, although DDSs are usually interpreted as extended 
dynamical systems with hyper-chaotic behavior, they ex- 
hibit distinctive features, including multiscaling and in- 
termittent driving fluctuations, not shared with other 
space-time chaos systems. 

A broad class of systems with delayed feedback are 
given by differential-delay equations such as 

y = ^ = ^{y^yr), (i) 

where yr — y{t ~ t) is the delayed variable and r is time 
delay. A canonical example is the Mackey-Glass model 
(sl (see below), for which it has been shown that the 
number of positive Lyapunov exponents grows linearly 
with the time delay ,17]. 

Generic surface growth model for DDSs.- Any DDS 
can be interpreted as a spatially extended system by in- 
troducing a simple transformation, t — x + 6 t, where 
X e [0, t] is th e sp ace variable, while S N is a discrete 
time variable |l8l ]. This is a powerful representation in 
which the formation and propagation of structures, de- 
fects and spatiotemporal intermittency can be more eas- 
ily identified 0] . Note that the time delay r becomes the 
system size in such a way that the time dependence with 
the delayed variable is transformed into a short-range in- 
teraction within the horizontal space coordinate x in the 
space-time representation. 

The dynamics of propagation of localized disturbances, 
i.e. the Lyapunov vector, can be obtained by linearizing 
(|H) which, after the transformation to space-time repre- 
sentation, becomes 19j 

y(j){x,9)^u(t){x,9)+v(j){x,9-l), (2) 
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where V = dx, 4>{Xi ^) = ^vi^^ ^) is the Lyapunov vector, 
and u = dT/dy and v = dT jdyT- 

Following Pikovsky et. al. Q, we can map this 
equation into a surface growth model by introducing the 
surface height as h{x, 9) — log \(p{x, 0)\. Then, after defin- 
ing dg(j){x, 9) = <j){x, 9 + 1) ~ (j){x, 9), one easily obtains 
the corresponding growth equation for the surface. We 
finally arrive to the equation. 



deh^~{l/v)Vh+{u/v) + l, 



(3) 



where V/i — d^h and the boundary conditions h{T,9) — 
h{0,9 + l). 

In the following we restrict ourselves to the most 
studied case of DDSs in the chaotic regime, in which 
u = const and w is a nonlinear function (see Eq.® be- 
low), so that 1/v will play the role of a noise source. 
As we shall see below, probability distribution and cor- 
relation of this noisy term will turn out to be essential 
for understanding the behavior of the growth equation. 
Note that, if Pv{v) is the probability density function 
(PDF) of V, the corresponding PDF oi ( = 1/v is given 
by P^(C) = |Cr^-Pt,(l/C)- Therefore, if P„ behaves as 
Pv{v) ~ Iv]"^ when u ^ 0, we have that the PDF of 
^ is then given by P(;(C) ~ for large |(| (con- 

versely for small \v\). This power-law decay distribution 
of the noise is very remarkable and leads to a dynamics 
dominated by rare fluctuations. The occurrence of large 
events of amplitude |C| ^ 1 has a small but non-negligible 
probability (they take place whenever the function v gets 
close to zero). As we will show below, these rare fluctu- 
ations drive the system towards a stationary state that 
displays multiscaling properties. Anticipating part of our 
numerical results, in Figure 1 we plot the PDF of the 
chaotic function ( for the Mackey-Glass model, which, 
as expected from our theoretical arguments, decays as a 
power-law P({C) ICI"^ ''^- 

The important fact that the surface growth equation 
that describes the propagation of disturbances in DDSs, 
Eq.Q, is dominated by rare fluctuations strongly sug- 
gests that the universality class of growth to which the 
surface must belong is the well-know linear Zhang model: 



dgh = DV^h + fo + Ti{x, 



(4) 



which is the representative model of the universality class 
of linear surface growth in the presence of power-law noise 
H EE El EE El' This model is expected to describe 
the dynamics of the Lyapunov vector in DDSs for large 
enough systems (i.e. long enough delay times). KPZ- 
like nonlinearities such as (V/i)^ are not expected to ap- 
pear in the effective surface growth equation, since Eq. 10) 
is statistically invariant under the h —f —h transforma- 
tion and this symmetry must be preserved. Moreover, 
for the class of systems we are interested in this Letter, 
u{x, 9) = const and a simple scaling argument shows that 
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FIG. 1: Numerical simulation of the Mackey-Glass model. In- 
set shows a typical realization of the chaotic function = 1/v 
in time, in which the pronounced non-Gaussian fluctuations 
are already visible. Main panel shows the power-law behav- 
ior of the probability density of obtained for time delays 
r = 500 (o) and r = 2000 (□). The line is a least-mean 
square fit to the data with slope —2.42. 



the additive noise in Eq.l^J must be ry oc v^- So, we ex- 
pect the effective noise term 77 in Eq. Q to be distributed 



according to a power-law, P{ri) 



where the 



noise index is fi = 2 {1 + v). Further justification that 
Eq.© is in the universality class of the linear Zhang 
model can be given as follows. By making use of standard 
stochastic functional techniques, like Novikov's theorem 
[20l[2l| . spatial average of the stochastic drift term (Vh 
can be calculated explicitely for the simpler case of Gaus- 
sian noise to obtain DV^h + bVh, where the constants 
D and b are related to space-time noise correlator inte- 
grals (see for instance Ref. l23| ) . In the case of noise with 
zero mean value and an even spatially homogeneous cor- 
relation function (77(2;, t)ri[x' , t')) one obtains a vanishing 
drift term 6 = 0. 

In the context of surface kinetic roughening, Eq.lffl 
corresponds to the well-known linear Zhang model [3, 
lla I14L lid |l6j that describes the surface height dynam- 
ics in the presence of an uncorrelated power-law dis- 
tributed noise source. In the following, we shall briefly 
review the scaling properties of this growth model. The 
height h(x^ t) exhibits scale-invariant properties char- 
acterized by the scaling behavior of the surface fluc- 
tuations as measured by the global width, W{L, t) = 

{[h{x,t)~h{^]^y/'^, where the average is calculated 
over all a; in a system of horizontal size L (overline) and 
noise (brackets). Space and time scale- invariance leads 
to the scaling behavior of the width. 



WiL,t) 



tP 



if i < 
if i > 



(5) 



The critical exponents a and P are the roughness 



3 



and growth exponent, respectively, and characterize the 
global scaling behavior of the surface height. The dy- 
namic exponent z is related with the horizontal extent 
of the correlation length ~ t^^'^ and fulfills the usual 
relation z = a/f3. A simple mean- field approximation 
0, Q provides estimates of the scaling exponents that 
depend on the index fi characterizing the power-law de- 
cay of the noise as follows. The dynamics is diffusive 
(z = 2) for any value of /i. However, the global rough- 
ness exponent depends on the noise index, a = (3 — 
for fj. < 2, when the fluctuations are dominated by 
rare events, and a = 1/2 for /i > 2, where the noise 
distribution tail becomes irrelevant and standard Gaus- 
sian behavior is recovered 0, Q|. Another interest- 
ing aspect of roughening dominated by power-law noise 
is multiscaling -namely, higher moments of the height- 
height correlation function do not scale with the same 
exponents- which indicates that the height distribution 
is far from Gaussian and dominated by rare fluctua- 
tions of the noise. In particular, multiscaling behavior 
in the linear Zhang model occurs for noise indexes in 
the range 1 < ^ < 3 123]. This behavior is analogous 
to that observed in the velocity fluctuations of turbu- 
lent fluids [23 as well as in some solid-on-solid models of 
epitaxial growth 0, 0] ■ Multiscaling properties can be 
proven by calculating the qth order height-height corre- 
lation function Gq{x,t) — {\h{xo + x,t) — h{xo,t)\'i)'''^'^ , 
which scales as Gq{x,t) ^ t^" for short times and satu- 
rates, Gq{x) ~ for long enough time 0,0,0] ■ 
Exponents aq and j3q that depend on the index q are the 
characteristic fingerprint of multiscaling behavior. 

Numerical results.- We have carried out extensive 
numerical simulations of several chaotic delayed systems 
in order to compare with our theoretical findings. We 
have studied the well-known class of systems 

y ^ -ay + F{yr), (6) 

which have been used in the past for a variety of appli- 
cations ranging from biology [8j to lasers with feedback 
I1II1II3- For the particular choices ^bp/{l+p^^), 
F{p) = b sm{p — po) and F{p) = b sm^{p — Po), which 
correspond to the Mackey-Glass Q, Ikeda and the 
delayed feedback model for secure optical communica- 
tions respectively. For this class of systems, the 
function u is constant, u — —a, and the additive noise in 
Eq.Q is ri{x,6) cx y/TJv, where the nonlinear function 
V = dF/dp. In all our numerical simulations we have 
used the Adams-Bashforth-Moulton predictor-corrector 
scheme |25j |. From the numerical solution of the corre- 
sponding equation for the evolution of the perturbations 
of the leading Lyapunov exponent, Eq.Q, we have con- 
structed the surface h[x,0) — \og\(j){x , 9)\ following the 
procedure above discussed to get the space-time map- 
ping. After that, the scaHng properties of the resulting 
surface have been studied. For the sake of brevity we 
focus the discussion of numerical results on the Mackey- 
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FIG. 2: Time evolution of the global surface width in the 
Mackey-Glass model for different system sizes L = 272, 5 x 
10^10^,2 X 10^,4 X 10^ (time delays). The straight line is 
a guide to the eye and has a slope (3 — 0.25. Inset shows a 
log- log plot of the stationary WaatiL) vs. L (fiat late regime 
of the main panel data). The line is a guide to the eye with 
slope a = 0.5. 

Glass model, but similar results were found for the other 
two systems studied. The parameters a — 0.1 and b — 0.2 
are used in all the results we are presenting here, and 
simulation with a time delay varying from ten to a few 
thousand time units have been carried out. The region of 
interest here corresponds to delays r 3> 17 for which the 
Mackey-Glass model is known to be hyper-chaotic [1 j . 

In Figure 1 (inset) we show a typical run of the sim- 
ulation for the Mackey-Glass system with a time delay 
T = 500 . The intermittent behavior of ( is already 
apparent in the form of large peaks in the inset plot. 
We recorded the values taken by v and collected statis- 
tics from long runs in order to obtain the PDF of the 
noise ^ — 1/v. In Fig. 1 we show our results for two 
instances of the time delay r — 500, 2000. We obtain 
a power-law distribution Pq{C) ~ ICI^^ "'^ that is inde- 
pendent of the time delay used and corresponds to an 
exponent v = 0.42 ± 0.02. The power spectrum of C (not 
shown) is flat indicating that ( lacks long-range tempo- 
ral correlations. We then expect the noise in the effec- 
tive surface growth model Q to be power-law distributed 
with ^ = 2 (1 + i^) = 2.84 ± 0.04. This value of the noise 
index p > 2 implies, according to the mean-field predic- 
tion for the linear Zhang model, that the roughness and 
growth exponent should be a = 1/2 and f3 = 1/4, respec- 
tively. Very interestingly, it also implies that one should 
observe multiscaling behavior, since p d (Ij 3]- 

In Figure 2 we plot the time behavior of the surface 
width for different system sizes L (i.e. time delays r) for 
the Mackey-Glass system. From that figure we obtain 
an estimation of the growth exponent f3 = 0.25 ± 0.03 in 
the intermediate times regime, before a stationary state 
is reached. The scaling regime appears after an earlier 
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FIG. 3: Multiscaling behavior of the stationary gth corre- 
lation function of the surface in the Mackey-Glass model 
(dashed lines) with a delay L = 10"^ for g = 1, 2, 3, 4, 5 (from 
bottom to top). Linear mean least-square fits (solid lines) 
give an estimation of the multiscaling roughness exponents as 
ai = 0.786(6), 02 = 0.496(2), 03 = 0.328(1), Q4 = 0.257(1), 
as = 0.210(1). 



The interpretation of DDSs as spatial systems has also 
been used in practical applications, in particular for the 
analysis of experimental data of CO2 lasers with optical 
feedback The technique allowed to identify many 

spacelike properties in experiments, including defects for- 
mation, transitions between weak and fully developed 
turbulence, etc. However, we have established a pre- 
cise difference between high-dimensional chaos in delayed 
systems and in truly extended systems (like coupled-map 
lattices, partial differential equations, etc). Our results 
also suggest that different high-dimensional chaotic sys- 
tems could be divided into a few universality classes 
based upon basic symmetries and conservation laws, akin 
to what occurs in scale-invariant surface growth. 
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from CONICET (Argentina) . Financial support from the 
MCyT (Spain) and FEDER under projects BFM2000- 
0628-C03-02, BFM2000-1108, BFM2001-0341-C02-02 as 
weU as project OCCULT IST-2000-29683 from EU are 
acknowledged. 



short transient regime (of about one decade), and be- 
comes longer for larger system sizes. For the largest 
system we used L = 4 x 10^ the scaling regime lasts 
for about two decades. After saturation, which occurs 
at times tsat ^ L^, the surface height becomes scale- 
invariant and the width grows as W{t tsat,L) ~ L", 
with a roughness exponent a = 0.50 ±0.05. From the in- 
set of Fig. 2, it becomes clear that our estimation of the 
roughness exponent fits better in the limit of large system 
size, showing that long delays are required to reach the 
true scaling regime. These exponent values agree with 
the theoretical values predicted for the Zhang model for 
a noise index \x = 2.84 > 2. 

Finally, we have also studied numerically the existence 
of multiscaling phenomena. In Figure 3, we plot the gth 
height-height correlation function for the Mackey-Glass 
system with time delay r = 10^. Multiscaling behavior 
is obtained for any time delay larger than r k, 17, which 
is the lower bound for chaos to appear in the system. 
However, multiscaling becomes clearer in larger systems 
(longer time delays), since the multiscaling region has 
a larger expand. As predicted by our theoretical argu- 
ments, multiscaling behavior is to be expected whenever 
II < 3. Crossover to self-affine scaling is observed on 
large scales as generically expected for surface roughen- 
ing driven by rare events (see 0, E| for details). 

Conclusions.- Symmetry and scaling considerations 
have led us to conclude that the linear Zhang model for 
surface roughening driven by power-law noise describes 
the dynamics of the Lyapunov vector for long enough 
time delays. We then argue that DDSs generically display 
multiscaling of fluctuations due to the power-law nature 
of the driving noise. 
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